Superconducting Coherence and the Helicity Modulus in Vortex Line Models 
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We show how commonly used models for vortex lines in three dimensional superconductors can be 
modified to include k = excitations. We construct a formula for the k = helicity modulus in 
terms of fluctuations in the projected area of vortex loops. This gives a convenient criterion for the 
presence of superconducting coherence. We also present Monte Carlo simulations of a continuum 
vortex line model for the melting of the Abrikosov vortex lattice in pure YBCO. 

PACS numbers: 74.60.-w (Type-II Sup.), 05.70.Fh (Phase Trans.), 75.40.Mg (Num. Simulations) 
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Phase transitions involving vortices in high temper- 
ature superconductors are the subject of intense study 
both experimentally and theoretically [0 . The enhanced 
thermal fluctuations strongly alter large parts of the 
mean field phase diagram [|2||3|], with new phases appear- 
ing, e.g., vortex line liquids, vortex glass phases, etc. A 
convenient quantity to study theoretically is the helic- 
ity modulus T (or spin- wave stiffness), which measures 
the free-energy increment associated with an externally 
imposed twist in the phase of the superconducting order 
parameter Q , and is proportional to the macroscopic su- 
pcrfluid density. Much work has been based on XY-\ikc 
models, defined in terms of this phase, where vortices 
appear only as topological defects. However, a formula- 
tion directly in terms of vortex degrees of freedom has 
several advantages. In this paper we show how the uni- 
form helicity modulus can be defined directly in terms of 
the vortex lines, without reference to the phase. We also 
show Monte Carlo results for T and other quantities in 
a continuum model of interacting vortex lines. 

One of the advantages of the vortex representation is 
the possibility to define the model on a continuum, and 
so avoid artificial pinning to a discretization lattice. Fur- 
thermore, one may include interactions coupling directly 
to the vortex lines, such as core energies and various 
types of disorder. The vortex representation therefore al- 
lows new parameter regimes to be reached compared to 
the phase representation. Both representations are fre- 
quently used in computer simulations |3]-[8) . T is straight- 
forward to calculate in the phase representation || . How- 
ever, in the vortex representation usually only k =/= fluc- 
tuations are included, making the uniform response T 
exactly zero, and extrapolations from finite k become 
necessary [g. Further more, when screening from gauge 
field fluctuations is taken into account, T(k) ~ k 2 for 
small, but finite k, since an imposed phase twist can be 
compensated by the gauge field. Extrapolation to k = 
then gives zero, but the small-k behavior of T(k) may 
still be related to the Meissner effect in a superconduc- 
tor and used to detect phase transitions [|| . The methods 
mentioned above involve extrapolations from the smallest 
available wave vectors to k = 0, thereby severely compli- 
cating the data analysis. An alternative may be to study 
winding number fluctuations, which are related to the 



magnetic permeability fx, but they suffer from being dif- 
ficult to equilibrate for large system sizes. 

In this paper we take a different route by modifying the 
vortex model in order to incorporate fluctuations with 
zero wave vector. The form of this modification is ob- 
tained using a duality transformation between the phase 
and vortex representations, paying due attention to the 
role of the boundary conditions. We show that periodic 
boundary conditions for the phases enter as an additional 
term in the Hamiltonian of the vortex representation. 
This allows direct evaluation of the k = helicity modu- 
lus in terms of fluctuations of the total net area of vortex 
loops, which can indeed be finite also in the presence of 
screening. The role of boundary conditions in the duality 
transformation has previously been explored in 2D lat- 
tice models f9Hli"| and 3D gauge glass models fL2f . Here 
we generalize this idea to continuous 3D systems, with fi- 
nite magnetic field, penetration depth, and temperature. 
Furthermore, we report on Monte Carlo simulations of a 
continuum London model. In contrast to previous con- 
tinuum simulations, which used 2D Bose models with 
planar interaction M, we take into account the full 3D 
long range interaction. Our model has the essential fea- 
tures to describe the vortex lattice melting transition in 
pure YBCO, where a continuum description should ap- 
ply- 

The starting point for our discussion is the Ginzburg- 
Landau theory in the London limit, where amplitude fluc- 
tuations of the superconducting order parameter 'F = 
|\&|exp(z0) are neglected. For simplicity we will use an 
isotropic continuum description, since our results are in- 
dependent of microscopic details. The generalization to 
other cases is straightforward. The Hamiltonian reads 
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where il — L x L y L z is the size of the system, 0(r) is the 
phase of the superconducting order parameter, B = Vx 
A is the magnetic flux density, H is an externally applied 
magnetic field, $o = hc/(2e) is the flux quantum, and 
J = <5q/ (167t 3 Aq) is a coupling constant with Ao the bare 
magnetic penetration depth. The first term in Eq. (jl]) 
is the kinetic energy with the superfluid velocity v = 
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V6* — (27r/$o)A, while the second describes the magnetic 
energy. The partition function is obtained by integrating 
over the phases 9(r), and gauge field A(r), subject to 
some gauge fixing condition: Z = J T>9T>' Ae~" H . In 
order to get a finite result a short distance regularization 
has to be imposed, e.g., by defining the model on a lattice. 
Physically this cutoff is of the order of the Ginzburg- 
Landau coherence length £o, and gives the size of the 
vortex cores. 

We now discuss the transformation of the model, 
Eq. (|l|), to a system of interacting vortex lines in some de- 
tail. For simplicity we start by considering the case with- 
out any externally applied magnetic field, H = 0. The in- 
teraction can be linearized by an integration over an aux- 
iliary field b(r), upon which the kinetic energy becomes 
J n d 3 rJ (ib ■ v + 5b 2 ). The superfluid velocity splits 
into a longitudinal part, describing the smooth spin- wave 
fluctuations of the order parameter, and a transverse 
part describing the singular vortices: v = vy + vj_. Inte- 
grating over the longitudinal part leads to the constraint 
Vb = 0, which can be enforced by setting b = V x a. Af- 
ter a partial integration of the first term and subsequent 
integration over B the Hamiltonian becomes 

H = jf d 3 rJ {27ml ■ m + 1 (V x a) 2 + ^ a^j , (2) 

where aj_ is the transverse part of a and m(r) denotes the 
vorticity, V x V9 = 27rm. Integrating out the auxiliary 
field a leaves us with the vortex Hamiltonian [l3|J^,p| 

H= I ^m(r) • V(r- r')m(r')d 3 rrf 3 r', (3) 

where K = (2ir) 2 J and V is the London interaction 



(4) 



k^ = 2irn fl /L fi , n M G Z, fj, = x,y,z. With periodic 
boundary conditions for the phases 9 we also get the con- 
straint of zero net vorticity J n m(r)<i 3 r = 0. 

In going from Eq. (|l]) to Eq. (|^) we implicitly as- 
sumed that v had no k = component, allowing us to 
throw away a surface integral. However, if the phases 
obey periodic boundary conditions, 9(r + L„) = 9(r) 
(where = L M e M ), there will be an additional energy 
term, coming from uniform fluctuations of v. An im- 
portant point here is that the integration over A should 
not include the uniform part Ao, since such fluctua- 
tions correspond to fluctuations in the boundary condi- 
tions. The additional energy is simply H' = ^ v 2 , with 
v = J a v(r)d 3 r. This is now related to the vortices as 
follows. The contribution to the vorticity m from a single 
vortex loop can be written 



m(r) = j> S(r - r')dr' = V x j 5(r - r')dS', 



(5) 



where T is a contour describing the vortex loop, and S 
denotes an oriented surface which has V as a boundary. 
Summing over all vortex loops gives the total vortex den- 
sity. Now, since the vorticity is the rotation of the super- 
fluid velocity, we may define 
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where the sum is over all vortices. This has the inter- 
pretation of the total projected net area of the vortex 
loops in each direction. Due to the periodic boundary 
conditions, the value of Q is uniquely determined by the 
positions of the vortices only up to an integer multiple 
of il/L^, reflecting the need to specify the total phase 
twist of the system. Thus, the variable Q keeps track 
of the total phase twist of the system and must be inde- 
pendently specified in addition to the vortices in order to 
completely specify the state of the system. The additional 
k = component of the energy is now given by 
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and the total Hamiltonian is given by the sum of Eqs. (||) 
and (|): H tot = H + H' . 

External magnetic field. — Assume now that n flux 
quanta <i>o penetrates the system in the z-direction. In 
this case the periodic boundary conditions for the phases 
have to be changed so that the system can accommodate 
a net number of vortices. For the vortices penetrating 
the whole system the area Q should now be measured 
with respect to a given reference line at some arbitrary 
but fixed position determined by the boundary condi- 



tions. A possible choice is to let 9(r + L p 



(r) 



rm — ^ J* r+ ** A • dr' , with the integral taken along a 
straight line across the system. In the gauge V • A = 
we may write A(r) = A per (r)-|- x (r — r ), where A per 
satisfy periodic boundary conditions and B = n^ li z /Q 
is the uniform part of the flux density. In this case the 
fixed reference line goes through ro in the direction of B. 

Helicity modulus. — The full importance of the new 
term in the energy becomes evident when one consid- 
ers the superfluid response of the system. Replacing 
the boundary conditions by twisted, #(r + L M ) — > 9(r + 
L M ) + 8, leads to the replacement — > — with 
= CIQ/L^ in Eq. (Q). This allows us to define the 
zero wave vector helicity modulus by 
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where ((Q 2 )) = (Q 2 ) - (Q) 2 , and F = -Tin Z is the free 
energy. T is non-zero in the superconducting state and 
vanishes at the phase transition. In the critical region 
of a continuous phase transition it obeys the Josephson 
scaling relation T ~ where £(T) is the correlation 

length a. 
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FIG. 1. Monte Carlo results for the (Ao,T) phase dia- 
gram for B = 0. The dotted line indicates the value of T c 
for the inverted Xy-model obtained in the limit Ao — * 0, 
T c ss Ka(Xo/a) 2 /3.0 (a is the lattice constant ~ £rj)- In the 
opposite limit, Ao — > oo, T c w 3.0A"a/(27r) 2 . Inset: T c is 
located at the intersection of curves for LT for different L. 
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FIG. 2. Helicity modulus T z and structure function S q at 
an ordering vector of the vortex lattice. Inset shows the jump 
in energy per vortex and layer at the transition. Also shown 
are two typical snapshots from the simulation, one below, and 
one above T c « 0.008Kd. 



The ordinary vortex Hamiltonian, Eq. (g), without the 
additional area term H' is recovered by integrating over 
and thus corresponds to having fluctuating boundary 
conditions 0. In this case the k = response is zero. 
Fluctuations in the winding number W = J n m(r)c? 3 r, 
are obtained only if the k = component of the magnetic 
flux density B is allowed to fluctuate. Then the mag- 
netic permeability fj, v = 4tt<9 (B v ) /8H v = ((B'l)) /flT = 
((W?)) §l/VLT can be calculated. 

Monte Carlo simulations. — To clearly demonstrate 
the practical usefulness of these ideas we now discuss 
Monte Carlo (MC) simulations. The inclusion of the area 
term in a simulation is straightforward. With the model 
defined in terms of vortices, the Monte Carlo moves con- 
sist of deformations of the vortex lines and (possibly) 
creations and destructions of closed loops. The change 
in the projected area coming from these local updates 
are accumulated in Q, and the change in total energy, 
including the area term Eq. (0), must be used to cal- 
culate the transition probabilities of the Markov chain. 
Optionally, these moves may be supplemented by global 
moves where Q M is changed by ifi/L^, corresponding to 
dragging a whole vortex across the entire system. The 
acceptance ratio for such moves can be expected to be 
quite low because of the high energies involved. Alterna- 
tively these global moves can be integrated out exactly, 
leading to the replacement of H' by a periodic Gaussian: 
exp(-/3ff') - Y.M^M-mrE^KiQ,- M,n/L^) 2 ) 
fl4| . Since this form of the area term leads to a somewhat 
more complicated expression for the helicity modulus it 
will not be used in what follows. 

We present simulation results for two different models: 
(i) a lattice superconductor in zero magnetic field with 
different values of Ao, (ii) a continuum vortex model of 



the melting of the Abrikosov vortex lattice. In both cases 
10 5 — 10 6 Monte Carlo sweeps were used, with the initial 
~ 10% discarded for equilibration. 

In Fig. |l| we present results from the first case, in which 
the phase transition is continuous. The critical temper- 
ature T c is determined from finite size scaling of the he- 
licity modulus Y(T, L) = L _1 T([T - T c ]L 1/ly ), as shown 
in the inset. Due to the new length scale given by Ao, 
scaling works only for rather large system sizes and cor- 
rections are clearly visible for small sizes, but the de- 
termination of T c is still quite accurate. In the limits 
Ao — > oo and Ao — > we recover known results, showing 
that our method works properly. Furthermore, a scaling 
collapse of MC data for A = 0.25 is obtained using the 
expected 3D XY value v m 2/3. 

In our continuum simulations in an applied magnetic 
field, we discretize the vortex lines only along the z- 
direction, using straight segments to interpolate between 
the zy-planes where the positions are continuous. We 
exclude overhangs and isolated loops, which should be of 
importance only close to the zero field T c . In addition 
to local MC moves of the positions in the xy-plane, we 
include moves where two flux lines are cut off and recon- 
nected to each other, allowing different permutations of 
the boundary conditions to be sampled. We use a full 3D 
long range interaction, given by Eq. (0) (with Ao = oo), 

supplemented by a Gaussian short distance cutoff e 4 ° , 
which acts between the midpoints of the vortex line seg- 
ments. The vortex lattice constant was set to 4£o and 
the layer separation d to 2£o- The number of layers for a 
system of N vortices was set to AyN. To avoid frustra- 
tion effects in the vortex lattice phase we use a hexagonal 
simulation cell with periodic boundary conditions in all 
directions. In Fig. we show the helicity modulus in the 
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direction of the applied field T z , and the structure func- 
tion Sq — (\m z (q)\ 2X j /(NL Z ) 2 at a reciprocal vector q 

of the Abrikosov vortex lattice, as a function of temper- 
ature. T x = T y = for all T, reflecting the absence of 
vortex pinning. At the transition both T 2 and Sq drop 
quite sharply, suggesting a first order melting transition 
to an entangled vortex liquid with no intermediate dis- 
entangled phase. Right at the transition, the time series 
of the internal energy or the structure function, obtained 
from the simulation, fluctuate around two different val- 
ues, giving further support for a first order transition. 
The inset shows how the average energy per vortex and 
layer approaches a jump at the transition as system size 
increases, with a latent heat of roughly 0.0015iVL 2 i ; ir. 
Taking into account the internal temperature dependence 
of the parameters [|15j, and using values for YBCO gives 
an entropy jump AS ~ 0.5/cb per vortex and layer, in 
rough agreement with experiments |l6| ] . 

We now comment on the implications of the ideas pre- 
sented above on the analogy between the statistical me- 
chanics of vortex systems and zero temperature quan- 
tum field theory of bosons in (2 + 1)1? The vortex 
lines are the analogue of boson world lines in a path in- 
tegral representation of the partition function, with the 
z-direction playing the role of imaginary time. A crystal 
ground state for the bosons corresponds to the Abrikosov 
vortex lattice phase in the vortex problem, while a su- 
perfluid boson ground state is mapped to an entangled 
vortex line liquid. The winding number fluctuations of 
world lines in (2 + 1)1? gives the superfluid density of the 
boson problem |L8| . It is of interest to study the conse- 
quences of our new area term H' in this context. To this 
end it is useful to view the London interaction, as be- 
ing mediated by a (massive in the screened case) (2 + l)D 
gauge field a M , see Eq. (||). The area term H' is now gen- 
erated by integrating over the k = component of the 
dual field strength / M = e^ vp d v a p (which corresponds to 
the auxiliary field b above Eq. (0) in the vortex problem) . 
By coupling the dual field strength, — (e v ,e x ,b z ), to 
an external source, <? M = (d y ,d x ,h z ), we find that the 
dual magnetic permeability ji = d (b z ) /dh z corresponds 
to the helicity modulus in z direction, T z |l]||. The in- 
verse dual dielectric constants e" 1 — d (ej) jddi in the x 
and y directions, correspond to the helicity modulus T M 
in the y and x directions, respectively. 

In summary, we have shown how to include k = fluc- 
tuations in vortex line models, and how the helicity mod- 
ulus can be obtained from fluctuations of the projected 
area of vortex loops. This is useful for detecting super- 
conducting phase transitions in models with and without 
screening of the London interaction. We also presented 
continuum Monte Carlo simulation results with a full 3D 
London interaction. This model should be appropriate 
for the vortex lattice melting in moderately anisotropic 
systems such as YBCO. Similar approaches should be 
useful in studies of, e.g., vortex glass transitions and 
quantum phase transitions. 
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